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Abstract 


This paper introduces a computational scheme for solving a class of aerodynamic 
design problems that can be posed as nonlinear equality constrained optimizations. The 
scheme treats the flow and design variables as independent variables, and solves the con- 
strained optimization problem via reduced Hessian successive quadratic programming. 

It updates the design and flow variables simultaneously at each iteration and allows flow 
variables to be infeasible before convergence. The solution of an adjoint flow equation 
is never needed. In addition, a range space basis is chosen so that in a certain sense 
the “cross term” ignored in reduced Hessian SQP methods is minimized. Numerical 
results for a nozzle design using the quasi-one- dimensional Euler equations show that 
this scheme is computationally efficient and robust. The computational cost of a typ- 
ical nozzle design is only a fraction more than that of the corresponding analysis flow 
calculation. Superlinear convergence is also observed, which agrees with the theoretical 
properties of this scheme. All optimal solutions are obtained by starting far away from 
the final solution. 

Key words, design optimization, constrained optimization, reduced Hessian meth- 
ods, quasi-Newton methods, successive quadratic programming 

AMS(MOS) subject classification. 65K05, 90C06, 90C30, 90C31, 90C90 

Abbreviated title. SCHEME FOR DESIGN OPTIMIZATION 

1 Introduction 

An aerodynamic design optimization problem can often be posed as 

min I(X,u) 

s.t. F(X,u) = 0, [ ’ 

where X £ denotes the discretized flow variables, and u £ 5? m denotes the design 
variables, which, for example, could be geometry parameters describing the shape of a 
profile; / : 3? n+m — ► ft is a cost function, which may, for example, measure the deviation 
from a desired surface pressure distribution; F : Ji n+m — ► $? n is a discretized version of the 
governing equations of the flow field. It is often the case that I and F are nonlinear, and 
the number of flow variables is much larger than the number of design variables (n >> m). 

Many computational methods for solving (1.1) have been developed, and each seeks to 
find the solution by solving a series of subproblems. One can categorize these methods by 
looking at how the design and flow variables are treated. In one category, design optimiza- 
tion approaches treat A as a dependent variable of u. Representatives in this category are 
the brute-force approach, the implicit gradient approach (e.g. [18]) and the adjoint equa- 
tions approach (e.g. [14]). These approaches have a feasibility requirement that is satisfied 
by solving the nonlinear flow equations to convergence in each design cycle. They only 
differ in how V U I is calculated at each iteration. The brute-force approach computes V U I 
by one-side finite differences, which requires m solutions of the flow equations. The implicit 
gradient approach and the adjoint equations approach also require at least two solutions of 
nonlinear flow equations at each design iteration. 
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In another category, design optimization approaches treat X and u as independent 
variables. Given X*, u*, each approach in this category finds AX*, Aw* and updates all 
the variables at once such that X* + AX*, u* + An* better solves (1.1), and iterates until 
a satisfactory solution is obtained. One interesting property is that the flow equations are 
not required to be satisfied until convergence. Approaches in this category are expected to 
mitigate the computational cost of repeatedly solving the nonlinear flow equations required 
by other methods. The success of flow solvers based on Newton’s method (e.g. see Barth 
[1] and [13]) also makes the all-at-once approach more tractable. 

The approach introduced in this paper belongs to this latter category. Methods that 
simultaneously update design and flow variables have been studied by many authors in- 
cluding Frank and Shubin [7], Huffman et al [13], and Hou et al [12]. The scheme in [13] is 
based on Successive Quadratic Programming (SQP) for equality constrained optimization. 
In their scheme the Hessian matrix is approximated by selectively calculating some second 
order terms instead of using an update formula. The calculation of second order terms 
could be prohibitive in many situations. The formulation introduced in [12] is based on 
successive linear programming. This scheme may suffer from convergence problems and 
in addition a linear adjoint system has to be solved at each iteration. Frank and Shubin 
[7] compare three optimization-based methods, including an all-at-once scheme (based on 
the full Hessian SQP that would be too expensive to form for large problems) for solving 
aerodynamic design problems. They point out, however, that an all-at-once scheme could 
be very efficient if carefully implemented. 

SQP is a mature and successful technique for solving nonlinear constrained optimization 
problems due to its relatively low computational cost and fast convergence. However, it is 
not widely used in aerodynamic design optimization. The major obstacle is that aerody- 
namic design optimization problems are often very large. In order to have an efficient and 
robust implementation of SQP methods, one has to take full advantage of special character- 
istics of the aerodynamic design problems. This paper is intended to address some aspects 
of this issue. 

Since the approach introduced in this paper originates from SQP for equality constrained 
optimization, we first describe the full Hessian SQP scheme in the context of (1.1). 

At each iteration of SQP methods, a quadratic approximation to the Lagrangian function 
of (1.1) 

L(X, u, A) = I(X , u) + A t F(X, u), 

is minimized subject to a linearization of the flow equations. This gives a subproblem 

min de»< n + m ) 9k d + \d T Hkd (1.2) 

subject to Fk + A^d = 0, (1.3) 

where d = ^ 9k = ^ ^ at ( Xk , Ufc), Ak = ^ ^ at (Xk,Uk), Fk = F(Xk, Uk ), 

and Hk is the Hessian (or approximation) of the Lagrangian function. 
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One difficulty implementing full SQP methods for aerodynamic design optimization is 
that the analytic Hessian of the Lagrangian is usually hard to obtain. Another difficulty 
is that approximating the Hessian by update schemes tends to create large dense matrices 
((n + m) x (n + m)). 

These difficulties can be overcome by dropping certain non-critical second order infor- 
mation. This idea has been pursued by many researchers including Coleman and Conn [4], 
Nocedal and Overton [16], Byrd and Nocedal [3]. We illustrate the basic idea as follows. 

One way of solving (1.2- 1.3) is via separation of variables. Suppose we can compute two 
matrices Yk E J£( n + m ) Xn and Zk E 9ft( n + m ) Xm such that the matrix [Yk Zk] is nonsingular 
and Z^Ak = 0 ( Zk is a basis for the null space of Aj). Let d = Ykd x + Z*d 2 and plug this 
into (1.2-1. 3). If we assume Z^HkZk is positive definite (which is reasonable because this 
matrix is indeed positive definite close to the solution due to the second order sufficient 
optimality condition), by standard techniques, d x and d 2 are given by 

<*i = -(AlY k )~ l F k 

d 2 = -{ZlH k Z k )-\tfig k + ZlH k Y k di). 

To avoid the computation of H k , the “cross” term Zj H k Y k di is simply ignored and 

H k Z k is approximated by an m x m matrix B k using a variable metric formula such as 
BFGS. Consequently, 


dx = -{A T k Y k )~'F k (1.4) 

d 2 = -B^ 1 Z k g k . (1.5) 

Methods based on (1.4- 1.5) are referred to as reduced Hessian SQP methods. 

The biggest advantage of reduced Hessian SQP methods for aerodynamic design opti- 
mization is that only a small m x m matrix holding second order information needs to be 
maintained and updated. Theoretically, reduced Hessian SQP methods are not as effective 
as full Hessian SQP methods. However, many versions are shown to have 2-step superlin- 
ear local convergence (for example see [16], [3] and [19]), which is good enough for most 
applications. 

One challenge of implementing (1.4- 1.5) for aerodynamic design optimization is the 
choice of Zk and Y*. Traditionally, Zk and Yk are chosen from an orthonormal matrix 
obtained from the QR factorization of A*. QR factorization based reduced Hessian SQP 
algorithms are proven to be very robust in practice (for example, see [11]). However, for 
aerodynamics design optimization, Ak contains the flow Jacobian matrix, which is usually 
very large. Orthonormal base are usually too expensive to compute. 

Many authors, including Gabay [8], Gilbert [9], Fletcher [6], and Xie and Byrd [19], have 
proposed more economical choices of Zk and Yk which are not orthogonal. One popular 
choice of Zk is 

Zt= (-(if>-W) .(-5), (1.6) 
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where, in case of aerodynamic design optimization, 5 is the sensitivity matrix. It can easily 
be verified that Zj Ak = 0. The calculation of the sensitivity matrix is usually affordable 
for many applications. An economical choice of Tit is [7 O] 7 ', which is used by many authors 
including Biegler et al [2], and Dennis et al [5]. 

This paper introduces a reduced Hessian SQP scheme for the aerodynamic design prob- 
lem (1.1). Equations (1.4- 1.5) are used to compute a trial step at each iteration. The 
popular choice of Zk given by (1.6) is used. However, a different choice of Y* is considered. 
We show that using [7 0] T has an undesirable effect of potentially producing a larger cross 
term. Here Yk is chosen such that Zj Yk = 0 holds, which in a certain sense will minimize 
the cross term. We also show that this choice of Yk has the advantage of gaining accuracy 
in the reduced Hessian approximation, which is very important for ensuring the reliability 
and robustness of a reduced Hessian SQP scheme. 

The new scheme is implemented and tested on the shape design of a transonic nozzle. 
The transonic flow through the nozzle with a given area ratio is governed by the quasi-one- 
dimensional Euler equations. The transonic conditions which generate a shock within the 
nozzle present a difficult test case, where methods typical of practical aerodynamic applica- 
tions are required. Such methods include, finite difference, finite element, and unstructured 
grid finite volume techniques employing various forms of highly nonlinear algorithm con- 
structions. 

This paper is organized as follows. First the aH-at-once reduced Hessian SQP scheme is 
introduced in Section 2. The solution of the quasi-one-dimensional Euler equations for flow 
through a nozzle with a given area ratio is discussed in Section 3. Section 4 first describes 
the design optimization problem used in our testing, then provides computational results 
and the performance evaluation of the new scheme on the design problem. The testing is 
carried out by using different numbers of design variables (spline coefficients for the nozzle 
shape) for the test problem. Finally, we give some concluding remarks in Section 5. 

2 An all-at-once reduced Hessian SQP scheme 

Let Bk be an approximate to ZjHkZk ■ Let Zk be given by (1.6). We define Yk by 



Clearly Z J Yk = 0 holds. 

Let d = Ykd\ + Zkd 2 . Replacing the reduced Hessian matrix by Bk and ignoring the 
cross term in (1.5), we have 

d 2 = -B^Z^gk- 

From (1.4), d x satisfies ( A[Y k )d x = -F k , which is equivalent to 
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which in turn is equivalent to 


J k (I + SS T )d, = - F k , (2.2) 

where «/* = (fx) (-Xfcj u k) is the current Jacobian matrix of the flow equation. Equation 

(2.2) can be solved by first solving J^y = -Fk for y, and then solving (I + SS T )di = y for 
d\. The solution of the former is simply the Newton step calculation from the flow equations 
at the current iteration, which can be provided by any Newton’s method based flow solver. 
The solution of the latter can be obtained by the conjugate gradient method, which is 
guaranteed to converge within (m + 1) iterations due to Rank(S S T ) = m (see Golub and 
Van Loan [10]). Another way of solving (I + SS T )di = y is by inverting (/ + SS T ) directly. 
It is easy to show that 

(/ + S5 T r 1 = / - 5(7 + S T S)- l S T . 

Note that (J + S T S ) is only an m x m matrix and its factorization can be obtained at 
minimal cost. 

After d\ and are available, d is given by 

d = Ykd\ + Zkd2 

= ( K s T ) d ' + ( / ) d2 ‘ 

Afterwards, we update our solution ^ ^ fc+1 ^ ^ + ad, update B k and the 

Langrange multiplier A* (which is used implicitly in the merit function calculation, see 
below), and go to next iteration. 

The Lagrange multiplier is asked to satisfy 

(Y k A k ) A* = -Y?g k , (2.3) 

which is equivalent to 

(/ + SS T )(J T \ k ) = -Y?g k . 

As a matter of fact, we will see that our scheme only needs the value of (J T A*). 

We establish the foDowing lemma that is useful to us. 

Lemma 2.1 Let M, N E $? nXm with n > m. If M and N have the same range space of 
dimension m, then 

M{N T M)~ l N T = M(M t M)~ 1 M t . (2.4) 
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(2.5) 


Proof. Clear ( N T M ) is invertible. Since M and N have the same range space, 

M(M t M)~ 1 M t = N(N t N)~ 1 N t . 


Postmultiplying each side by M(N T M) 1 N T and using (2.5) gives 

M(M T M)~ 1 M T M(N T M)~ 1 N T = N(N T N)~ 1 N T M(N T M)- 1 N T 
=» M{N t M)~ x N t = N{N t N)-'N t 


which gives (2.4) using (2.5). □ 

When certain update criteria are satisfied, the Hessian approximation update is carried 
out via the BFGS formula 


Bk+i = Bk - 


BkSkS^Rk ykvl 
slB k Sk y^sk' 


( 2 . 6 ) 


where 


Vk = Z% +l [V (x,u)L(Xk+i,Uk+i, \k) - V(J )c,u)L{Xk,Uk,\k)] (2.7) 

$k = (^k+ i^fc+i) (2-8) 


with a being the step length resulted from a line search. This choice of y k , $k pair is 
among the choices recommended by Nocedal and Overton [16]. From the definition of the 
Lagrangian function 


Z J+\[{9k+\ + -Afc+lA*;) — ( 9k + -Afc-Xfc)] 


z k+il9k+i - ( 9k ~ A k (Y? Ak)~ x Y l fc T < 7 *)] 

(2.9) 

Zk+il9k+i - (I - Y k (Y? Yk)-'Y?)gk] 

(2.10) 

Z k+\[9k+\ - Z ki Z k Z k)~ l z k 9k\- 

(2.11) 


Equation (2.9) to (2.10) is due to Lemma 2.1. Equation (2.10) to (2.11) is because 

Z k {ZlZ k )-'Zl + Yk(Y k T Yk)- 1 Y? = I 


due to Z^Yk = 0. 

To ensure convergence, a merit function is needed to monitor the progress towards the 
solution. We choose to use the /i merit function for its simplicity and low computational 
cost. The l\ merit function is defined as 

The directional derivative of 4>^ along d is given by 

D K k (x > u ; d ) = sl d - rII^IIi- (2.12) 
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Plugging d = Ykdi + Zkd 2 into (2.12) yields 

D<t,^ k (X,u;d) = glY k d x + Z k d 2 - /*fc||/fc||i. 

Clearly, on the one hand 

9I ^kd 2 = —g£ ZkB k l Zk gk < 0, (2.13) 

since Bk is forced to be positive definite. On the other hand, from (2.3) and d\ = 

-(Y?Ak)-'F k , 

glYkd , = A lF k = {J T \k) T {J~ l Fk). (2.14) 

Note that both (J T Xk) and ( J~ l Fk ) are available from previous calculations. From (2.13- 
2.14) 


D <t>» k (X,u;d) < 0 

can be guaranteed by choosing pk > (A^ F* : )/||Fjt|| 1 

We point out that J T is not needed in the entire scheme, which is desirable in aerody- 
namic calculations. 

One feature of our choice of Yk is that the cross term ignored in the reduced Hessian 
SQP formulation is likely to be minimized. Since this is not a obvious claim, We establish 
the following lemma. 

Lemma 2.2 For a given matrix N G & nXm and a vector p G 9? m , over all the choices of 
matrices M G JJ nXm such that (N T M) is invertible, \\M(N T M)~ 1 p\\ is minimized if M has 
the same range space as N . 

Proof. Let N = U DV T be the singular value decomposition of N with V G U nXm and 
D y V G SR mXm , and M = QR be the QR decomposition of M with Q G 5? nxm and R G 
jjmxm. Clearly, D and R have full rank. Then 

||M(1V T M)- I p|| = \\QR(VDU t QR)-'p\\ = \\Q(U t Q)- 1 D- 1 V t p\\ 

= \\{U T Q)- l D-'V T p\\ > ||D -1 V rT p||/CT mai , (2.15) 

where a max is the largest singular value of ( U T Q ). Since U and Q are unitary, ||17 T Q|| < 
||tf r ||||£2|| = 1, which implies a max < 1. Hence 

||M(^V T M)- 1 p|| > \\D-'v T p\\. 

The equality is achieved when U T Q is unitary (in this case <r max = 1. Hence it is sufficient 
to show that U T Q is unitary if M and N have the same range space. 

Clearly if M and N have the same range space, then Q and U have the same range 
space. Since Q and U are unitary, QQ T = UU T holds. Therefore 

(U t Q) t (U t Q) = Q t (UU t )Q = Q t (QQ t )Q = ( Q t Q)(Q t Q ) = /, 
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which implies that U T Q is unitary. □ 

Recall that d\ = — {A^Yk)~ l Fk- Choosing Yk such that Z^Yk = 0 implies that Yk and 
Ak have the same range space due to Zj Ak = 0. Hence by Lemma 2.2, ||yjt<fi|| achieves 
minimum value when Z k Yk = 0. For a given Zk and Hk, a smaller ||Fjfcdi|| is desirable since 
it has a major contribution to the magnitude of the cross term Z^HkYkd\ ignored in the 
reduce Hessian formulation. 

Choosing Yk such that Z^Yk = 0 also has the potential advantage of improving accuracy 
in the Hessian approximation. Note that (2.8) implies 

Zk+\$k = Zk-\-\{Z'k +x Zk+\) 1 Zk + \ad. 

The Taylor expansion of yk used in the step secant update gives 

yk = Zk +l V(x ,u)(X,u)L(Xk, Uk, \k)<xd + 0(||ad|| 2 ) 

= Zl +i ^lx,u)(X,u)L(Xk, Uk, Xk)Zk+i{Zk + i Zk+i ) _1 Zj +1 ad 

+Zl+i^fx,u)(x,u)L(Xk,‘Uk, Xk)(I - Zk+i(Zk +1 Zk+ 1 ) 1 Z][ +1 )ad 

+o( IMI 2 ) 

= Zk+\ ^(X,u)(X,u)L(X k , Uk, Xk)Zk+lSk 

+ Zk+ iV ( 2 Xt u){X,u)L{Xk, Uk, (Yk+\Yk +\ ) l Yk+\<xd + 0(||ad|| 2 ). 

Since Bk+i is expected to approximate the reduced Hessian matrix, it is desirable to keep 
the size of the second term as small as possible. This goal can be partially achieved by 
controlling the size of Vit+i (Yf +l Yfc+i ) _1 Y^ d. Using d = Ykd\ + Zkd 2 , we have 

iKn+i^in+o-'n+iMii 

= miYfYkr'Yfw+omi 2 ) 

= \\Y k dt + (Yk(Y k T Yk)- l Y k T )Z k d 2 1| + 0(||d|| 2 ). (2.16) 

Choosing Y k such that Zj Y k — 0 makes || || achieve its minimum value as analyzed 

above, and makes ||(F* ; (y fc T Yk ) ~ 1 Y k )Zkd 2 \\ zero. Hence it is likely that this choice would 
yield the lowest value for (2.16). 


3 Quasi-one-dimensional Euler equations 

For our purposes here, we have chosen one popular form of central finite differences with 
nonlinear artificial dissipation applied to the quasi- ID Euler equations. 

The quasi- ID Euler equations are 


where 


F(Q) = d x E(Q) + H(Q) + D(Q) = 0 0.0 < x < 1.0 



pu I 0 

pu 2 + p , H = -pd x a(x) 

_u(e + p)J 0 


(3.1) 


P 

pu 

€ 


E — a(x) 
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(3.2) 



with p (density), u (velocity), e (energy), p = (7 - l)(e - 0.5 pu 2 ) (pressure), 7 = 1.4 (ratio 
of specific heats), and a(x ) = (1. - 4.(1 - at)x(l - x)) (the nozzle area ratio), with a t = 0.8 
. For a given area ratio and shock location (here x = 0.7) an exact solution can be obtained 
from the method of characteristics. 

We choose one popular form of central finite differences to discretize these equations. 

d x q « S x qj = - j = 1 , . . . , 7 max 

Ax = 1.0 /(Jmax — 1), Uj = u(jAx) (3.3) 

It is common practice and well known that artificial dissipation must be added to the 
discrete central difference approximations in the absence of any other dissipative mechanism, 
especially for transonic flows, see Pulliam[17]. 

For simplicity here, we use a constant coefficient dissipation of the form 


D 4 ( Q) = V x (6< 4) A x V x A r Q J ) 

with 


(3-4) 


~ 9j- 1 * = q j+1 - q 3 (3.5) 

with a typical value of = Ag. 

Boundary condition at j = 1 and j — J max are defined in terms of physical conditions 
(taken from exact solution values) and will be treated as Dirchlet (fixed conditions) for now. 
The total system we shall solve is 


The Jacobian matrix for the linear 


dissipation case of Eq.(5) can be written exactly as 



[C] 

[D] 

[E] 



[B] 

[C] 

[B] 

[E] 


[A] 

[B] 

[C] 

[D] 

A - d -L- 

dQ- 



[Ah 

Wi 


v 


\ 


[E] 






[Ch 

[D)j 

[E], 




[A] 

[ B] 

[C] 

[D] 

[E] 



[A) 

[5] 

[C] 

[D] 

[E] 



[A] 

[B] 

[C] 

[D] 




[A] 

[B] 

[C] / 
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where the elements ([/] is the 3 X 3 identity matrix) are defined as 

Mt = 

IC\, = + {CJ\, 

ID), = -4<W[/] + 

[E\, = <* 4 >[/] 

with 

0 1 


[AJ]j = o(*)i 


0 

(7 — 3)u 2 /2 —( 7 — 3)« (7-1) 

[ [-7 u e/p + (7 - l)ti 3 ] [7«/P~ f (7 — 1)« 2 ] 7« 


(3.8) 


and 


[CJ]j = -d r o(x)j 


0 0 0 

(7-l)tr 2 /2 -(7-l)« (7-1) 

0 0 


(3.9) 




The d x a(x)j term is done with central differences and there are some slight modifications 
to the dissipation terms near the boundaries. 

4 Test results for a nozzle design problem 

The nozzle design problem. We assume that a target velocity distribution u*, is given 
for each computational grid. The design problem we are trying to solve is 

Find j/i, i = 1, • • - ,m (spline coefficients describing a(x)), such that 
2 12j=i X ( u j ~ u j ) 2 ts min imized subject to (3.6) being satisfied. 

For our test examples, the breakpoints of the spline are evenly distributed in the interval 

[ 0 , 1 ]. 

Some implementation issues. Each design is started from a flat nozzle and a flat flow 
solution. The first reduced Hessian approximation is given by B 0 = Zq Z 0 . 

A trial step d is calculated as described in Section 2. A quadratic backtracking line search 
is carried out using the /1 merit function, which yields a step length a*- The following update 
criterion is used for the reduced Hessian approximation. At iteration fc, is updated using 
the BFGS formula (2.6) if 


SkVk > 0. lotfcll Vfcrfi |[. 


(4.1) 
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Otherwise, we set Bk+\ = Bk . This criterion is devised to guard against loss in accuracy 
of 5*, as recommended by Xie and Byrd [19], and is critical for ensuring convergence of a 
reduced Hessian SQP algorithm. 

We give test results on design problems with zero, two and ten design variables. Through- 
out this section the number of grids Jmax is set consistently to 67. We chose this number 
since having an even breakpoint distribution for the splines requires that ( Jmax — 1) be 
divisible by the number of design variables plus one. Hence each problem has 67 x 3 = 201 
flow variables, and the same number of constraint equations. The convergence was moni- 
tored by ||JV^|| + ||F fc || (||F*|| in case of zero design variables), where Nk = Zk(Z][ Zk)~ l Z^ 
is a projector onto the null space. This quantity indicates how well the first order necessary 
optimality conditions are satisfied, and should be zero at the true solution. The tolerance 
is set to 10 -7 . 



Figure 1: Computational results without design variables, (a) Target (dash line) and final 
(checked solid) flow solutions, (b) Convergence of norm of residual versus iterations. 

Test results without design variables. These results are obtained by calling the flow 
solver with the optimal nozzle area ratio. We include these results as a baseline for mea- 
suring the cost of our optimization scheme. These results are also useful for examining the 
effect of the presence of design variables on the flow solution. The overshooting in the final 
solution can be controlled by nonlinear dissipation, as described in Pulliam [17]. 

Figure 1 (a) gives the final and target flow solutions. Figure 1 (b) shows the convergence 
history of the nonlinear flow residual versus the iterations. Figure 2 gives snapshots of partial 
solutions at various iterations. It can be seen that the transition is quite violent due to the 
existence of a shock in the solution. 

Test results with two design variables. For this problem four spline coefficients are 
used. We keep the coefficients at the end points fixed, leaving the area ratios at the two 
between points as design variables. Figure 3 gives the test results for this case. Figure 3 (a) 
shows the initial, final and target flow solutions. Figure 3 (b) shows the initial, final and 
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Figure 2: Snapshots of flow solutions at various iterations without design variables Dashed 
line - target solution; Checked solid line - current solution. 


target nozzle solutions. Figure 3 (c) gives the convergence history of the objective function. 
Figure 3 (d) plots the convergence history of the quantities ||./Vjk<7fc|| and ||F^||. 

The optimization converged quite rapidly taking 14 iterations. Figure 3 (a) and (b) show 
that the final flow and nozzle solutions match well to the corresponding target solutions. 
Figure 3 (c) shows that interates converged to a Kuhn-Tucker point (point that satisfies 
first order necessary conditions) superlinearly. 

Test results with ten design variables. For this problem twelve spline coefficients are 
used. Again we keep the coefficients at the end points kept fixed, leaving the area ratios at 
the ten between points as design variables. Figure 4 gives the test results for this case. 

It took 30 iterations for the optimization process to converge to the given tolerance. 
Figure 4 (a) and (b) show that the final flow and nozzle solutions match well to the corre- 
sponding target solutions. 

The convergence is slower than the case with two design variables. We believe that this 
is partially due to over-parameterization that makes the design problem much harder to 
solve. Nevertheless, as shown in Figure 4 (d), it converged rapidly (superlinearly) to the 
solution after wandering around for a relatively longer transient period. 

For the situation with two design variables, it was observed that flow solutions went 
through a less violent transient regime compared to the situation when running the flow 
solver without the presence of design parameters. For the situation with ten design variables, 
similar phenomena occurred. Figures 5 and 6 give some snapshots of intermediate flow 
and nozzle solutions with ten design variables at various intermediate iterations. These 
snapshots indicate that the transient is much milder with the presence of design variables. 
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Figure 3: Computational results with 2 design variables, (a) Target (solid line), initial 
(dashed line) and final (checked solid) flow solutions, (b) Target (solid line), initial (dashed 
line) and final (checked solid line) nozzle profiles, (c) Reductions in objective function, (d) 
Reduction in projected gradient (+) and flow residual (*). 


We suspect that this is partially due to the fact that an intermediate partial nozzle solution 
actually helps the flow solver to locate the shock and to estimate its magnitude. 

Cost comparisons. Finally, we discuss the performance of the reduced Hessian SQP 
scheme with respect to the pure analysis (solution of nonlinear flow equations). Table 1 
indicates that optimization with two design variables is only slightly more expensive than 
the analysis (costing 1.37 nonlinear flow solves). For the same problem, a brute force finite 
difference approach would require 3 nonlinear flow solves for a single optimization iteration. 
To solve the design problem, the number of nonlinear flow solves required by either an 
implicit gradient method or an adjoint equations approach would be at least twice as much 
as the number of optimization iterations which may easily exceed 10. For the situation with 
ten design variables, a fairly tough problem due to over-parameterization, the performance 
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Figure 4: Computational results with 10 design variables, (a) Target (solid), initial (dashed) 
and final (checked solid) flow solutions, (b) Target (solid), initial (dashed) and final (checked 
solid) nozzle profiles, (c) Reductions in objective function, (d) Reduction in projected 
gradient (+) and flow residual (*). 

of the reduced Hessian SQP scheme is still reasonably good. It costs an equivalent of 
8.217 nonlinear flow solves, which would compare favorably over the other three methods 
mentioned above for similar reasons. 

5 Concluding remarks 

An efficient and robust reduced Hessian SQP scheme for aerodynamic design optimization 
has been introduced in this paper. This scheme requires neither the solution of an adjoint 
equation nor any exact second order information. The optimization was able to converge 
from tough starting points. Computational results show that this scheme has a potential 
to achieve the same order of cost as one nonlinear flow solve. However, much more work 
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Figure 5: Snapshots of flow solutions at various iterations. Solid line - target solution; 
Checked solid line - current solution. 


Number of Design Variables 

Total Flops 

Number of Iterations 

Cost Ratio 

0 

2142176 

20 

1.000 

2 

2935839 

14 

1.370 

10 

17602144 

30 

8.217 


Table 1: Efficiency of the all-at-once scheme 


is still needed before we have a practical scheme. We believe the efficiency of this scheme 
can be improved if combined with grid sequencing [20] or multigrid techniques [15]. These 
techniques give efficient ways of providing a good starting point for the design on the final 
mesh. Issues such as dealing with inequality constraints and applications to 2D and 3D 
problems are currently under investigation. 
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